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We report the study of a model of a two-level system interacting in a nondiagonal way with a 
complex environment described by Gaussian orthogonal random matrices (GORM). The effect of 
the interaction on the total spectrum and its consequences on the dynamics of the two-level system 
is analyzed. We show the existence of a critical value of the interaction, depending on the mean level 
spacing of the environment, above which the dynamics is self- averaging and closely obey a master 
equation for the time evolution of the observables of the two-level system. Analytic results are also 
obtained in the strong coupling regimes. We finally study the equilibrium values of the two-level 
system population and show under which condition it thermalizes to the environment temperature. 
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I. INTRODUCTION 

The nonequilibrium statistical mechanics of small quantum systems has become a topic of fundamental importance 
for nanosciences. It is indeed very important to understand what the minimum sizes and conditions arc under which a 
quantum system can display a relaxation to an equilibrium state. Isolated and finite quantum systems have a discrete 
energy spectrum, which has for consequence that all the observables present almost periodic recurrences on long time 
scales. Nevertheless, their early time evolution may still present relaxation types of behavior that is important to 
study. Although tools have been developed in nonequilibrium statistical mechanics to describe such relaxations to a 
state of equilibrium by master and kinetic equations, the conditions of validity of these equations remain little known. 

It is the purpose of this paper to contribute to the clarification of these questions of validity of the kinetic description 
by studying of a simple model of a two-level system or spin coupled to a complex environment described by random 
matrices. Indeed, work done during the last decade has shown that the Hamiltonian of typical quantum systems 
presents properties of random matrices on their small energy scales. Here, we consider the Hamiltonian of the 
environment as well as the operator of coupling between the spin and the environment to be given by random 
matrices taken in a statistical ensemble of Gaussian orthogonal random matrices. This defines a model in which many 
results can be obtained analytically. 

Similar models using Gaussian orthogonal random matrices 0, Q or using banded random matrices 0, 0, IH 
have been studied. In |5(, it has been shown that random matrices used as environment coupling operators have a 
universal feature. 

The model we here consider differs from the spin-boson model by the density of states of the environment. Instead 
of monotonously increasing with energy as in the spin-boson model, the density of states of the environment obeys 
Wigner's semicircle law in our model and is thus limited to an interval of energy with a maximum density in between. 
We notice that such densities of states appear in systems where a spin is coupled to other (possibly dissimilar) spins 
as, for instance, in NMR, in which case the density of states of the other spins forming the environment also present a 
maximum instead of a monotonous increase with energy. Our model may therefore constitute a simplification of such 
kinds of interacting spin systems. Our main purpose is to understand the conditions under which a kinetic description 
can be used in order to understand the relaxation of the spin under the effect of the coupling with the rest of the 
system, which we refer to as a complex environment. 

The plan of the paper is the following. The model is presented in Sec. [H] The properties of the spectrum are 
described in Sec. IIIII The relaxation in the time evolution of the spin is studied in Sec. |IV| The very long time 
behavior and the approach to the equilibrium is discussed in Sec. [V] Conclusions are drawn in Sec. IVII 

II. THE MODEL 

We are interested in the study of a total system composed of a simple system (with a few discrete levels) interacting 
with a complex environment (with many levels). Wc consider a two-level system as a prototype for the simple system. 
For this kind of total system, the time-dependent Schrodinger equation is of the following type: 
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where 

• o'k) o'j/) and o"2 are the 2x2 Pauli matrices, 



<j z is the Hamiltonian of the two-level system, 



• 

• A is the energy spacing between the two levels of the system, 

• Hb is the Hamiltonian of the environment, 

• a x is the coupling operator of the system, 

• B is the coupling operator of the environment, 

• A is the coupling parameter between the system and the environment. 



The well-known spin-boson model j^l lJ llfl | is a particular case of the total system where Hb corresponds to an 

infinite harmonic oscillator lattice and B is linear in the degree of freedom of the environment. 

Here, we want to define a new model, the spin-GORM model, also described by the Hamiltonian Q and for which 

Hb and B are Gaussian orthogonal random matrices (GORM) (see appendix A for some basic property on GORM). 
Let us discuss now the spin-GORM model in more detail. 

As we said, we want to model a two-level system that interacts with an environment that has a complex dynamics. 
Here, complex is used in a generic way. The complexity can come, for example, from the fact that the cor resp onding 
classical system is chaotic like in a quantum billiard or for the hydrogen atom in a strong mag netic field [lllll2j |. It can 
also come from large coupling in an interacting many-body s yste m as in nuclear physics [ill or in interacting fermion 
systems such as quantum computers [TH IT^ | . Wigner in 1960 14, 15, 16] was the first to develop random-matrix theory 
for the purpose of modeling spectral fluctuations of complex quantum systems containing many states interacting with 
each other. This tool has now become very common in many fields from nuclear physics to quantum chaos. This is 
the reason why we consider random matrices to characterize the complexity of the environment operators. 

The environment operators of the spin-GORM model, Hb and B, are defined by 

Hb — °~ndX 

B = a§ D x', (2) 

where X and X are two different y x T Gaussian orthogonal random matrices with mean zero. Their nondiagonal 

(diagonal) elements have standard deviation a^D = 1 ( a D = V^)- X an( l X are two different realizations of the 

same random matrix ensemble and have therefore the same statistical properties, c^f) and cr^ D are the standard 

deviations of the nondiagonal elements of Hb and B, respectively. For these random matrices, the width of their 
averaged smoothed density of state is given by 

VH B = cr% B n V8N, 



' ND 



VB = a% D V8N (3) 
(see Appendix A). 

It is interesting to define the model in such a way that, when N is increased, the averaged smoothed density of state 
of the environment increases without changing its width T>Hb- The width can be fixed to unity. This is equivalent 
to fixing the characteristic time scale of the environment. For doing this, it is necessary to rescale the parameters as 
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follows: 



a = a% B D - 
t = at , 
A = ^ 



'8N 



A = A^ , 



N = N. 



The time-dependent Schrodinger equation of the spin-GORM model becomes 

: d|*(t)) 



ih- 



dt 



(4) 



with the rescaled total Hamiltonian 



= H s ~ 
A 

= 2°' 



H B + Xa x B 
1 



=A + A(j T 



/8iV 



(5) 



As announced, we have now T>Hb = T>B = 1. 

In the following, without loss of generality, a will always be taken equal to unity. Notice that, to model an 
environment with a quasicontinuous spectrum, the random matrix must be very large (TV — > oo). 

In order to get ensemble averaged results, one has to perform averages over the different results obtained for each 
realization of Eq. 10. When we use finite ensemble averages, the number of members of the ensemble average will 
be denoted by x- 

We see that the Hamiltonian JSJ) is characterized by three different parameters: A, A, and N. We define three 
different parameter domains in the reduced parameter space corresponding to a fixed iV in order to facilitate the 
following discussion. These three regimes are represented in Fig. ^ domain A with 1 > A, A; domain B with 
A > 1, A; domain C with A > 1, A. 

III. THE SPECTRUM 

In this section we study the spectrum of the complete system for the different values of the parameters. This study 
is important in order to understand the different dynamical behaviors that we encounter in the model. 

Let us begin defining the notations in the simple case where there is no coupling between the two parts of the total 
system (A — > 0). The isolated system has two levels separated by the energy A: 



Hs\s) = sj\s), 

where s = ±1. The environment has the standard spectrum of a GORM 

6g\b)=E?\b), 



(6) 



(7) 



where b — 1,2, ...,N/2. The Hamiltonian of the total system without interaction between the system and the envi- 
ronment is thus 



Hn — Hs 



H F 



(8) 



and the spectrum is therefore given by 



H \n) 



E°\n) 



(9) 
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with n = 1, 2, TV and 



The eigenvectors are tensorial products of both the system and environment eigenvectors: 

\n) = \s)®\b). (11) 

Let us now define the notations in the opposite simple situation where the coupling term is so large that the 
Hamiltonian of the system and of the environment can both be neglected (A — > oo). Using the unitary matrix U 
acting only on the system degrees of freedom 



U 



l l 

\/2 V2 



the total Hamiltonian becomes 

= Xa z B. (12) 
E KTI and \nrj) = \k) ® \rf) are, respectively, the eigenvalues and eigenvectors of the Hamiltonian 

Ho\kti) = XE K7I | K-q) = XkE v \kti), (13) 

where 77 = 1, ^ and k = ±1. 

After having defined the notation in the two extreme cases A — > and A — * 00, we will start the study of the 
spectrum with interaction A 7^ 0. 

The total spectrum is given by the eigenvalues {E a } that are solutions of the eigenvalue problem: 

H tot \a) =E a \a), (14) 

where a = 1,2, ...,N. It is very difficult to obtain analytical results for this problem. We will therefore study the 
total spectra using a method of numerical diagonalization of the total Hamiltonian. 

A. Smoothed density of states 

In order to have a quantitative understanding of the global aspect of the spectrum (on large energy scales), we will 
study the total perturbed averaged smoothed density of states. 

The environment-averaged smoothed density of states obeys the semicircle Wigner law [see Eq. ()82|l in Appendix 
A]: 



AN / 1 1 
n-(e) = V V ( 2 )2 - £2 if |£|< 2 

= if |e| > ~, (15) 

where e is the continuous variable corresponding to the environment energy E^ . 

Therefore, when A = 0, the total averaged smoothed density of states is the sum of the two environment semicircle 
densities of states which correspond to both states of the two-level system [see Eq. JSJ]: 

n(e) = n w (e-j)+n w (s + j) 



AN [TZ " A~7 AN [TZ " A~7 

- VV ( 2 )2 - (£ -2 )2 + VV ( 2 )2 - (£+ 2 )2 



lf (a" 2 )<|6|<( 2 + ^ 

elsewhere (16) 
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where e is the continuous variable corresponding to the total energy E®. The semicircle densities n w (e — and 
n w (e + y) are schematically depicted in Fig. [21 for different values of A. The numerical density of states of the total 
system corresponding to A = is depicted in Fig. [3] 

When A — > oo (meaning that the coupling term becomes dominant in the Hamiltonian) , the averaged smoothed 
density of states of the total system (see Ea. l|13|) ') is given by 

„(£)=„»(£)+„«(_£) = ^(^) 2 -(e) 2 if |e|<^ 

= if |e| > ~ (17) 

where e is the continuous variable corresponding to the total energy E KJ1 . This result can be observed in Figl^a) for 
A = 10 (because A, 1 < A). 

When A ^ 0, the total averaged smoothed density of states is also plotted in Fig. |3| The main observation is that 
there is a broadening of the complete spectrum when one increases A. In Fig. |3Ja), we see that the averaged smoothed 
density of states changes in a smooth way from l|16fl to i|17fl . But in Fig. OJb) and much more in Figs. Ijjc) and (d), 
the two semicircle densities n w (e — — ) and n w (e + — ) seem to repel each other as A increases. This is due to the 
fact that the levels of a given semicircle density don't interact with each other but only interact with the levels of the 
other semicircle density. This is a consequence of the nondiagonal form of the coupling. Therefore, having in mind 
the perturbative expression of the energies [see Eq. (|95|l of Appendix C], one understands that when A is nonzero, 
the eigenvalues that are repelling each other with the most efficiency are the ones closest to the center of the total 
spectrum. 



B. Eigenvalue diagrams 

The global effect of the increase of A on the eigenvalues has been studied with the average smoothed density of 
states. But in order to have an idea of what happens on a finer energy scale inside the spectrum, it is interesting to 
individually follow each eigenvalue E a as a function of A on an eigenvalue diagram. 

The first thing to note (see Fig. is that the increase of the coupling induces a repulsion between the eigenvalues. 
This has, as a result the broadening of the total spectrum as we already noticed on the smoothed averaged density 
of states. If one looks closer inside the fine structure of the eigenvalue spectrum, we see that there is no crossing 
between the eigenvalues. This is a consequence of the fact that there is no symmetry in the total system. Therefore, 
the non-crossing rule is always working. Each time two eigenvalues come close to each other, they repel each other 
and create an avoided crossing. One can notice that there is a large number of avoided crossings inside the region 
where the two semicircle n w (s — — ) and n w (s + — ) overlap (see Fig. [3 in order to visualize the overlapping zone 
that extends from — — 1 to 1 — — ). But outside this overlapping zone, there are very weak avoided crossings (and 
of course no crossing) and all the eigenvalues appear to move in a regular and smooth way. The regions seen in 
Figs. Ela), (b) and (d) are inside the overlapping zone and the one in Fig. 0Je) is outside. In Fig. 0Jc), the lower 
part of the diagram is inside the overlapping zone and the upper part is outside. This phenomenon is due to the 
fact that the eigenvalues of a given semicircle density do not interact with the eigenvalues of their own semicircle 
density but only with those of the other semicircle density. It is the consequence of the nondiagonal nature of the 
coupling in the spin degrees of freedom. When the coupling becomes large enough (A > 1), the avoided crossings also 
disappear inside the overlapping zone. Around A = 1 and inside the overlapping zone, there is a smooth transition 
from an avoided-crossing regime that gives rise to a turbulent and complex A evolution to another regime without 
much interaction between the eigenvalues that gives rise to a smooth A evolution. 



C. Spacing distribution 

The eigenvalue diagrams only give us a qualitative understanding of the fine energy structure of the spectrum. For 
a more quantitative study, it is interesting to look at the spacing distribution of the spectrum. 

When A = 0, each semicircle distribution, corresponding to a different system level, has a Wignerian level spacing 
distribution 

p-( s ) = | se -f* 2 . (18) 

But the total spectrum is made by the superposition, with a shift A, of two of such semicircle distributions. The 
shift creates a Poissonian component to the total spacing distribution in the overlapping zone {% — 1 to -f + 1). A 
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Poissonian distribution is given by 



P*{ 3 ) 



= e 



— s 



(19) 



Therefore, we choose to fit the total spacing distribution by the mixture 

pf u (s) = dP w (s) + (1 - Ci)F p (s). 



(20) 



This choice of the form of the fit is empirical but reasonable because the correlation coefficient of the fit is always 
close to one (between 0.963 and 0.982). 

We computed the spacing distribution and made the fit 1)20(1 in order to compute the mixing coefficient C\ for 
different values of A in the case where the overlapping zone covers almost the whole spectrum (A <C 1). The results 
are plotted in Fig. \5\ We see that there is a specific region of A values where the total spacing distribution is close to 
a pure Wignerian one. This region corresponds to the situation where X 2 N = Oil), i.e., when the typical intensity of 
the interaction between the nonperturbed levels is 0(A 2 ) (since the first nonzero correction in perturbation theory is 
of second order due to the nondiagonal coupling in our model), becomes of the order of the mean level spacing O(jj) 
in the total system. In this region, the level repulsion is maximal and effective among almost all the states. 

In the limit A — > oo the spectrum is again the superposition of two semicircle distributions, one for the Ei^s and 
one for the E-i n 's according to Eq. (fO|l . As a consequence, one gets a Poissonian type of spacing distribution. 

The other A cases will have a Wignerian spacing distribution outside the overlapping zone and a mixed one (like 
for the case A -C 1) inside. This can be visually seen on the eigenvalue diagram that we studied before. 



We want now to have some information about the eigenstates of the total system inside the overlapping zone of the 
two semicircle densities. Following [(|, we define the quantity 



If we fix e° and study e°) as a function of e, we will call it the local density of states (LDOS). 

If we fix e and study £(e, e°) as a function of e°, we will call it the shape of the eigenstates (SOE). 

We here focus on the SOE. The SOE tells us how close a perturbed eigenstate (A ^ 0) at energy e is from the 
nonperturbed eigenstate (A = 0) at energy e°. If the SOE is a very narrow function centered around e = e°, the 
concept of a nonperturbed eigenstate is still useful. This regime corresponds to very small coupling, for which the 
interaction intensity is lower than the mean level spacing < X 2 N < 1, and will be called to the localized regime. 
In the limit N — > oo this regime disappears. If one increases the coupling, the interaction intensity between the 
nonperturbed states begins to be larger than the mean level spacing between the states: A 2 iV > 1. The nonperturbed 
levels start to be "mixed" by the interaction and the SOE starts then to have a Lorentzian shape with a finite width T, 
centered around e = e . This regime is called the Lorentzian regime. If one further increases the coupling parameter, 
the SOE begins to spread over almost the whole spectrum. This regime is called the delocalized regime. 

In the banded random matrix model of 6] , the regimes are classified according to a different terminology and there 
is an additional regime corresponding to a spreading that goes beyond the energy range where the coupling acts (due 
to the finite coupling range of the banded matrices) . The motivation for our change of terminology will become clear 
in the study of the dynamics. 

The Lorentzian regime can be separated into two parts. For small coupling, the width of the Lorentzian T is smaller 
or of the order of magnitude of the typical energy scale of variation of the averaged smoothed density of state of the 
environment 8e: n(e + Se) w n(e). For larger coupling, the Lorentzian width extends on an energy scale larger than 
the typical energy scale of variation of the density of states of the environment. Therefore, one has r < Se in the 
former case and T > Se in the latter case. The different regimes are represented in Figs. El and 

To represent the SOE of the spin-GORM model, we discretize the energy axis in small cells of the order of the mean 
level spacing and we average the SOE over x realizations of the random matrix ensemble. We see in Fig. [HJa) the 
typical shape of the SOE going from the perturbative regime (A = 0.01, 0.05) to the beginning of the Lorentzian one 
(A = 0.1). In Fig. IHIb), we see the SOE across the Lorentzian regime (0.2 < A < 0.8). We also see the delocalized 
regime, when A — > oo and the SOE gets completely flat (A = 10). Figure[8Jc) shows the width of the Lorentzian from 
a fit made on the SOE curve. The correlation coefficient of the fit helps us to determine the region of the Lorentzian 
regime where the SOE is very well fitted by a Lorentzian. It has been verified that the width of the Lorentzian 
is independent of N in the Lorentzian regime. Figure Efd) (log-log) shows that the width of the Lorentzian has a 
power-law dependence in the coupling parameter close to two in the Lorentzian regime. 



D. The shape of the eigenstates (SOE) 




(21) 



a , n 
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E. Asymptotic transition probability kernel (ATPK) 

An interesting quantity, which is close to the SOE, but which has a nice physical interpretation, is the asymptotic 
transition probability kernel (ATPK) . 

The transition probability kernel (TPK) gives the probability at time t to be in the level |n) if starting from p(0). 
It is defined as 

n t (n|p(0)) = (n\e- i6Mt p(0)e i6Mt \n) (22) 
The ATPK is the time average of the TPK 

OooM/KO)) - km - f dtU t (n\p(0)) 

= ^|(a|n)| 2 (a|p(0)|a). (23) 

a 

The distribution of the ATPK in energy is given by 

iMeVo)) = Y, noo(n|p(0))£(£„ - e°). (24) 

n 

The ATPK has an intuitive physical interpretation. It represents the probability after a very long time to end up in 
a nonperturbed state |n), having started from the initial condition p(0). The ATPK is the convolution of the SOE 

n oo (e°|e ') = ££Ka|n>| 2 Ka|m>| 2 a(£ n - e a )5{E m - e '), (25) 

n,m a 

with p(0) = £ m \m)(m\S(E m - e '). 

Because we are interested in a random matrix model and for the purpose of studying the dynamics, we will perform 
averages of two different kinds. The first kind of average is a microcanonical average over states belonging to the same 
given energy shell of width Se for a given realization of the random matrix ensemble used in our total Hamiltonian. 
The second kind of average is an ensemble average over the x different realizations of the random matrices ensemble. 
For the microcanonical average, a choice of the width of the energy shell Se has to be done in such a way that it is 
large enough to contains many levels (to get a good statistics) and small enough to be smaller than or equal to the 
typical energy scale of variation of the averaged smoothed density of states of the environment Se. Therefore, the 
adequate choice of the width of the energy shell corresponds to -h < Se < Se. 

Different ATPK are depicted in Fig. (SJ where there is no random matrix ensemble average x = 1- I n Fig- EI a ); the 
ATPK is a microcanonical average inside the energy shell at energy e° and of width Se° . "Pin" denotes the total 
probability of staying inside the energy shell after a very long time. A A 2 = 1 is on the border between the localized 
and the Lorentzian regime and A = 1, 10 in the delocalized regime. We see that, until AA 2 = 1, the main probability 
stays in the initial energy shell. But when the Lorentzian regime starts, the probability spreads over energies larger 
than Se (here Se » Se° ). In Figs. GJb) and (c), no average has been done. The initial condition is a nonperturbed 
pure state (corresponding to an energy close to zero) and one can see the individual probability of being on another 
nonperturbed state for the different regimes. We can see that in the localized regime the probability of staying on 
the initial state is much more important than the probability of leaving it. We also see that the Lorentzian regime 
starts when the initial state loses its privileged position containing the main probability and, therefore, when the 
neighboring levels begin to have an important fraction of the total probability. 

IV. TIME EVOLUTION 

We now want to understand the time evolution of our model and more specifically the population dynamics of the 
system. 

The exact evolution of the total system is described by the von Neumann equation 

i>(t) = -i[H tot ,p(t)], (26) 

where iftot is given by Eq. (J5J. The system dynamics is obtained from (|2(j|l using the reduced density matrix 
ps(t) = Tr Bp{t). The total system has a finite and constant energy. At initial time, the environment has a given 
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fixed energy e corresponding to a microcanonical distribution inside a energy shell centered at e and of width Se. Se 
is chosen in such a way that it is large enough to contains many levels (to get a good statistics) and small enough 
to be smaller than or equal to the typical energy scale of variation of the averaged smoothed density of state of the 
environment Se. Therefore, the adequate choice for the width of the energy shell corresponds to < Se < Se. The 
dynamics is then averaged over the x realizations of the random matrix ensemble. In the following, we consider in 
detail the two different extreme cases of weak and strong couplings. 



A. The weak coupling regime (A <C 1) 

We derived in Ref. |l7j a perturbative equation for the description of the evolution of a system (with a discrete 
spectrum) weakly interacting with its environment (with a q uasicontinious spectrum). This equation has been shown 
to be equivalent to the well-known Redfield equation |10L Il8t Il9| when the typical energy scale of the system (typical 
energy spacing between the system levels) can be considered small compared to the typical environment energy scale 
Se (typical energy scales on which the smoothed density of states of the environment varies). In the spin-GORM 
model, this condition means that n(e + A) sa n(e). 

We will not perform in this paper the detailed derivation of this equation and its application to the spin-GOE 
model. This has been done in Ref. (l7|. We will simply recall the main ideas and results of this paper and apply 
them to the study of the dynamics of our model. 

The main idea is to suppose that the total density matrix can be described at all time by a density matrix of the 
following form: 

p(t) = -^—J2\ s ^ s '\ p ^(ti B ;t), (27) 

n(H B ) , 



where the environment density of states is 



n(e) = Tr B 6(e - H B ), (28) 



Doing this, we neglect the contribution to the dynamics coming from the coherence of the environment but we keep 
those of the system. Therefore, the total density matrix is diagonal in the environment degrees of freedom but not in 
the system ones. For comparison, let us recall that in the derivation of the well known Pauli equation, both types of 
coherences are neglected and the total density matrix is completely diagonal j^jj, El • 
The matrix of elements P ss '(H B ;t) is Hcrmitian, 

P ss ,(H B ;t) = P;, s (H B ;t). (29) 
The reduced density matrix of the system becomes 

Ps(t) = Tv B p(t) = J deTv B S(e-H B )p{t) = \ s )( s '\ J deP ss ,(e;t). (30) 

s.s' 

We see that each element of the reduced density matrix of the system depends on the environment energy. This is 
fundamental in order to take into account the finite energy effects of the total system. 

Using Eq. I|27|l and performing a perturbative expansion up to the second order in A on Eq. (|26[) (see Ref. |17|1. 
one gets for the population dynamics the equation 



P ss (e;t) = -2A 2 ^ [ ( S \S\s')(s'\S\s)P Ss (e;t) J de'F(e,e')n(e') 



, sm(E s - E s 



E, - Eg, +e-e> 



-(s\S\s)(s'\S\s)n(e) f de'F(e,e')Pgg,(e';t) S[li{ ^ jf g ' + 6 f N , (31) 

where F(e,e') = "| (e|_B|e') | 2 " , where the quotes denote a smoothening over a dense spectrum of eigenvalues around 
e and e'. We see that the probability Pg g (e;t), if initially concentrated in a given energy shell, can spread in energy 
under the dynamics. This is a typical non-Markovian effect happening on a short-time scale. It is due to the presence 
in the equation of the energy integral and of the sm ^ T - ) function that has a finite width in energy at short time. This 
equation is non-Markovian in the sense that the coefficients of the differential evolution equation are time-dependent 
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and that this time dependence can be neglected on long-time scales (performing the Markovian approximation) when 
the environment has a faster dynamics (i.e., typically a larger energy scale) than the system. 

The Markovian approximation consists of taking the infinite-time limit of the time dependent coefficients using 
the property lim T ^ 00 sln ^ T ) = tt6(£). This approximation is justified if the contribution of the £ = value (if it 
exists) has the main and almost unique contribution to the energy integral. If one further neglects the contributions 
of the coherence to the population evolution (this is automatically satisfied for the spin-GORM model because of the 
nondiagonal coupling), one gets a Pauli-type equation 



P ss (e;t) = -2ir\ 2 £ \(s\S\s')\ 2 F(e, E s - E s , + e)n(E s - E s , + e )P ss (e;t) 



s'^s 



-2tt\ 2 \{s\S\s')\ 2 F(e,E s - E s , + e)n(e)P s , s ,(E s - E s , +e;t). 

s'^s 



(32) 



We see that the Markovian approximation strictly keeps the dynamics of the total system inside an energy shell. 
Starting with the probability located on a given energy shell, the dynamics will preserve the probability inside this 
shell. But of course, the probability of the different states inside the shell are varying. The dynamics described by 
this equation can be seen as a random walk between nonperturbed states of the total system belonging to the same 
energy shell with transition rates between these states given by the Fermi golden rule. 

We now apply our equation to the spin-GORM model in order to study of the population evolution through a z 
(the difference between the probability of being in the upper state of the system minus the probability of being in the 
lower one). Doing this, we suppose that the environment is quasi-continuous N — > oo and that the random matrix 
ensemble average has been performed x ~ * 00 ■ F° r the non-Mar kovian equation (|31|l . one gets 



\NM 



(*) 



de' [P ++ (e';t)-P__(e'+A;t)]. 



where 



\ 2 r + 5 

P±±{e]t) = — I de 

7T 



\ z . ,sin(±A + e - e')t 



(±A + e - e<) 

In the Markovian limit, Eq. (|33|) becomes 

(a z ) M (t) = P ++ (e;t)-P^(e + A;t), 
which obeys a Pauli-type equation so that we get 

(&*) M (t) = (&*)£+[(&,) M (0) e-* 
where the equilibrium value of the populations is given by 



\-{e) 2 -J\-{e + AY 



\- (e) 2 + J\-(e + AY 



and the relaxation rate by 




(33) 

(34) 

(35) 
(36) 

(37) 
(38) 



where we adopted the convection that \fx — if x < 0. 

In the case of the spin-GORM model, the transition probability between states belonging to the same total energy 
shell only depends on the smoothed density of states. This is due to the fact that the environment coupling elements 
between the environment nonperturbed states are randomly distributed because the environment coupling operator 
is a random matrix. This has the consequence that the equilibrium value of the populations l|37|l corresponds to a 
microcanonical distribution probability of the states belonging to the total energy shell independently of the initial 
distribution of these states inside this energy shell. 

We notice that, in the general case, the Markovian assumption made on our equation does not directly give a Pauli- 
type equation which is an equation for the populations only. To get a Pauli equation, the further approximation, 
which consists of neglecting the contribution of the coherences to the population dynamics, has to be done. For the 
particular case of the spin-GORM model, this further approximation is not necessary because it is automatically 
satisfied. 
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B. The strong coupling regime (A ^> 1) 



We are now interested in describing the dynamical regime where the coupling parameter A is very large and, 
therefore, dominant in front of 1 and A. We will again suppose that the environment is continuous N — > oo and that 
the random matrix ensemble average has been performed \ ~ > 00 ■ 

The population dynamics of the system is given by 



(a z )(t) = Tr p{t)o z 
Using the following unitary transformation acting on the spin degrees of freedom 



U 



l l 
75 "75 



we get 

(a z )(t) = Tr e l{ ^^+" B+x ^ 6)t U^ps{Q)UpB^)e- l{ ^^ +&B+x ^ &)t a x 
Using the following perturbative expansion of the evolution operator to order zero in 4: 



Or, 



we get 

(*«>(*) 

Using the following notation 
and 



Tr e iXt& ' B Wps{0)Up B (0)e 



i\ta z B 



O 



a z B\nrj} = a z \K,)B\rj) = kE v \kt]} 



H B \b)=E B \b), 



we find that 



(a z )(t) A =° ^e iXtKE ^{ V \p B (0)\v){K\Wps(0)U\-K)e iXtKE -+O 



K , 1 1 



Because ps is diagonal in the basis that diagonalizes Hb-, we have 

MM0)to>=£|<#>| 2 <&|pi»|&>. 



(39) 



(40) 



(41) 
(42) 

(43) 

(44) 
(45) 

(46) 
(47) 



If we perform an ensemble average over different realizations of Hb and use the random matrix eigenvectors statistics 



it we pcrioi 



we find 



(v\pb(o)\v) (Hb) = £|(#>| 2 ^ B >|/M&> 

b 

V- 2 (Hb) 2 



(Hb) 



-(Hb) 



N 



Equation 146|) therefore becomes 
( Hb ) i->o 2 



<**>(*) 



N 
2 

N 



J2e 2tXtKE v (k\U^ Ps {0)U\ - k) 

,n 

(11^^(0)^71 - 1) ^e 2 ^ + (-l\mp s (0)U\l))J2< 



(48) 



(49) 



-2i\tE„ 
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Performing now the following ensemble average over different realizations of B: 

-(B) 



,2i\tE„ 



de n[e)e 



2i\U 



AN 



de J--e 2 e 2lXU 
nli V 4 



= N 



Ji(Xt) 

xt ' 



(50) 



we finally get 



(**>(*) 



(H B ,B) i^O J! (At) 

At 



<l|t>-+p s (0)t>| - 1) + (-11^(0)^11) 



(51) 



It is easy to show that the term of order j is zero. This explains the fact that O (^) has been replaced by O (p-)- 



Choosing as an initial condition 



we get 



Mo) 



1 
0/' 



(52) 



(53) 



We have found a well-defined behavior of the system dynamics when the coupling parameter A is so large that the 
coupling term can be considered to contribute alone to the dynamics. 



C. Numerical results 



We will now numerically study the validity of the approximated equation that we just derived in order to understand 
the dynamical evolution of our model. 

In our numerical simulations, the initial condition of the system is always the upper state (|52|) : (<x 2 )(0) = 1 . The 
different parameter domains represented in Fig. ^wiU play a fundamental role in our discussion of the dynamics. 

We begin by the comparison between the results of the non-Markovian version (|33|l and (|34|l of our perturbative 
equation and the results of the Markovian version H36(l - (|38|l (of Pauli type) in order to understand better the conse- 
quences of the Markovian approximation. Figures ITUT a') and llOf b) show the time evolution of (ct z ) for both equations 
at different values of the coupling parameter. The time axis has been scaled by the coupling parameter (A 2 i). This 
scaling is characteristic of the Lorentzian SOE regime. The values of A have of course to be reasonably small in 
order to remain consistent with the fact that these equations are obtained perturbatively. The time scale that we are 
observing is the global one: from the initial condition to the equilibrium. The characteristic energy of the system (A) 
is different in Figs. HOt a) and llOf bl. But we are in both cases in domain A of the reduced parameter phase space (see 
Fig. nj. For such small values of A, we see almost no difference between Figs. EJ a ) and llOf bl. The characteristic 
time scale of the environment is of order 2tt, and is therefore in both cases much shorter than the system one. We 
are in situations where the Markovian approximation makes sense on times longer than 2ir. Using the X 2 t scaling, 
the Markovian equation is independent of A. This is not the case for the non-Markovian equation. We see that the 
stronger A is, the larger the deviation is between both equations. This is the consequence of the fact that the X 2 t 
scaling that we use amplifies the non-Markovian short time behavior (that occurs on time of order of the characteristic 
time scale of the environment 2ir) when A increases. Figs. HOf c) and ll(Jf d) show us the long time behavior of the 
Markovian and non-Markovian versions of our equation. We see that the equilibrium value of (a z ) depends on A for 
the non-Markovian equation. This is not the case for the Markovian equation. The differences between the equilib- 
rium values are small. But, comparing Fig. llUf with llOf dl. one notices that when A is larger, the differences are 
larger. This is a consequence of the error made on short-time dynamics using the Markovian approximation. Because 
this error is more important when A is large, the consequence on the long time dynamics is more important, even if 
globally small. 
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We now compare the non-Markovian version (|33fl - l|34[) of our perturbative equation with the exact von Neumann 
equation (|26[) . staying in the parameter domain A. Figures ITTT a) and II If show the time evolution described by 
both equations at different A values using the X 2 t time scaling. In Fig. Illf aV the Markovian and the non-Markovian 
versions of our equation are so close that we only plotted the second one. We see that the non-Markovian equation 
is valid not only for values of A below an upper bound, but also above a lower bound. When A is too small, the 
non-Markovian version of our perturbative equation does not fit with the exact result. The exact initial dynamics is 
well reproduced by our equation, but the relaxation process to the equilibrium value is not reproduced. This is due 
to the discrete nature of the spectrum and therefore depends on the number N of states and disappears in the limit 
N — > oo. It corresponds to the border between the localized and Lorentzian regimes of the SOE. This is one of the 
main results of this paper. This phenomena is of course again related to the effect of the perturbation between the 
levels in the total spectrum that we already observed in the spacing distribution, in the SOE, and in the ATPK. A has 
to be large enough (A 2 > -h) to "mix" the nonperturbed levels in order for all the states inside the total unperturbed 
energy shell to be mixed together. Remember that the equilibrium value of our Markovian perturbative equation is 
given by a microcanonical distribution in the total unperturbed energy shell (see Eq.J2ZJl). One also sees in Fig. Illf b) 
that our perturbative equation loses again its validity above a certain value of A. It happens at value of A that cannot 
be considered as perturbative any more. It also corresponds to a value of the coupling parameter corresponding in 
the SOE to the transition in the Lorentzian regime when the width of the Lorentzian starts to be larger than the 
typical variation energy scale of the environment density of states. 

In Fig. Illf c'). we deal with the parameter domain B (see Fig. We again compare the non-Markovian version 
of our perturbative equation with the exact dynamics given by the von Neumann equation. In this case, the system 
dynamics (of period OAtt) is faster than the environment one (of period 2ir). We are in a highly non-Markovian 
situation. The Markovian version of our equation (that describes no evolution in this case) completely misses the 
observed behavior of damped oscillations. But the non-Markovian equation reproduces this behavior with a very 
high accuracy. The frequency of the oscillations corresponds to the system dynamics ones and the damping of these 
oscillations occur on a time scale corresponding to the environment characteristic time scale. 

Finally, in Fig. Illf d). we are in the high coupling parameter domain C (see Fig. QJ. We see that when A becomes 
large enough to make the coupling term dominant in the total Hamiltonian, the dynamics obey the Bessel behavior 
derived in Ea. (|53() . It is important to notice that this behavior scales in time according to Xt. 

A summary of the validity of the different approximated equations (Markovian (|36[) - i|38|) and non-Markovian l|33|) 
and 134|) versions of our perturbative equation and strong coupling Bessel equation H53|) ~) and of the different scalings, 
depending on the regime that one considers, is represented in Fig. 1121 

D. Average versus individual realizations 

An interesting point is the comparison, for the system dynamics, between the averaged (random matrix ensemble 
averaged or microcanonically averaged) dynamics and the dynamics of an individual realization within the statistical 
ensembles. This latter corresponds to a dynamics generated by an initial condition that corresponds to a pure state 
and without any random-matrix average (x — !)■ 

We see in Fig. EH c ) _ (f) the dynamics of a system with small energy spacing A = 0.1 for different values of A 2 iV. 
The solid line represents the random-matrix and microcanonically averaged dynamics. The dashed lines depict some 
of the individual members of the random-matrix ensemble corresponding to an initial pure state inside the total 
unperturbed energy shell. We see that the larger A 2 iV is, the closer the individual trajectories are from the averaged 
trajectories. We therefore have, when A 2 iV is large enough, that the individual realizations are self-averaging in 
the random-matrix and microcanonical ensembles. In order to quantify this behavior, we plotted in Fig. 113( a) the 
variance between the individual trajectories and the averaged trajectories as a function of time for different values of 
A 2 TV. We observe that this variance decreases as A 2 TV — > oo. In Fig. I13f . we show that the asymptotic value of the 
variance decreases with a power-law dependence with respect to A 2 TV. 

We again see the relation with the SOE regimes. In the localized regime (Fig. I13f c)). the dynamics is governed by 
very different individual trajectories oscillating with a very few frequencies that differ from one individual trajectory to 
another. This is a consequence of the fact that the perturbed levels are still close to the nonperturbed ones and are only 
slightly affected by neighboring nonperturbed levels. In the Lorentzian regime (Figs. I13f d) and (e)), each individual 
trajectory follows roughly the averaged trajectory and contains a very large number of different frequencies. This 
shows that the interaction "mixes" many of the nonperturbed levels, deleting the discrete structure of the spectrum. 

Therefore, we can say that our master equation (|31fl or (|32|l holds with a given accuracy for a majority of individual 
trajectories if A is small enough satisfying A > CN~ V with v < | and a constant C > 0, in the limit TV — > oo. 
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V. THE VERY LONG TIME BEHAVIOR 

We here focus on the very long time behavior of our model, in other words on its equilibrium properties. 

A. Equilibrium values of the system observables 

Let us consider the spin observable a z of the two-level system. This observable evolves in time according to 

(a z )(t) = Tr p{Q)e lflMt a z e- tflMt 

= Y,^W)W)Wz\a)e i ^- E ^ t . (54) 

ex., a.' 

The time-averaged value of (a z )(t) is obtained performing the following time average: 

i <- T 



(o- z )oo = lim - / dt{a z )(t) 

= £<a|p(0)|a><a|S->>. (55) 

a 

We see that this time averaged value clearly depends on the initial condition. 

An important result, that holds everywhere on the parameter space and for all kinds of initial conditions, is that 
the observed equilibrium value given by the exact von Neumann equation corresponds very well to the time-averaged 
value (|55|l . Therefore, the study of the equilibrium properties of systems as ours is equivalent to studying the time 
averaged quantities 

Let us note that the standard initial condition we used till now is a microcanonical distribution around energy e 
for the environment and an upper state for the system formally given by 

P(0) = HXH = £ ^^M<H (56) 



Therefore, we have that 



(a|p(0)|a) = £ S{E " £) |Hln)| 2 . (57) 

n ^ ' 

An interesting point is to understand when the Markovian perturbative equation gives the correct equilibrium value. 
We want, therefore, to compare (|55|l with 1)3 7|) . This is done in Fig. I14lwhere we plotted the time averaged value (<7z)oo 
as function of the initial energy of the environment e (the initial value of the total system being l|5t)[) ). The different 
curves in Figs. fMT a') and IT4T b) correspond to different values of A. We compare these curves to the A independent 
curves given by the Markovian perturbative equation l|37|) . As expected from our precedent study of the dynamics, 
we find again that, when A is too small, the "mixing" between the nonperturbed levels is not sufficient inside the 
microcanonical energy shell and the Markovian perturbative results overestimate the equilibrium values. If A is too 
large, the Markovian perturbative equation gives again bad results. The Markovian perturbative equilibrium values 
are correct in a characteristic region of A. The beginning of this region corresponds, in the SOE, to the critical value 
of A at which the transition occurs from the localized to the Lorentzian regimes. The end of this region corresponds 
to A values that cannot be considered any more as perturbative. Figures H4T c') and !14f d) show that (p z )oo scales like 
A 2 A. This again confirms our precedent analysis. 



B. Thermalization of the system 

One of the important questions is to understand the conditions under which the system thermalizes under the effect 
of a weak contact with the environment or in other words, under which conditions the system relaxes to a canonical 
distribution corresponding to the microcanonical temperature of the environment. 

We begin by recalling these conditions in the general case of a small system weakly interacting with its environment. 
The isolated total system is composed of the system and the environment and has the total energy 



#tot = e + e. 



(58) 
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e is the energy of the environment and e the energy of the system. We suppose that the contribution of the interaction 
energy between the environment and the system is negligible compared to the total energy. The microcanonical 
environment entropy is defined as 

S B (e) = klnSl B (e), (59) 

where f2s(e) is the number of states of the environment available at energy e. This number can be related to the 
density of states of the environment n B (e) using the fact that £l B (e) = n B (e)Se, where 5e is a small energy interval 
but contains many states of the environment. The microcanonical temperature of the environment is given by 

1 _ dS B (e) {m) 



T B {e) de ' 

It can be expanded in the system energy as 

T B (e = E tot -e) = T B (e = E tot ) - = Etot) + . . . , (61) 

because we suppose that the system energy is much smaller than the environment energy. The specific heat capacity 
of the environment is 

1 (62) 



C vB (e) de 
If the condition 

\C vB (e)\ » 



(63) 



T B (e) 

is satisfied, the temperature expansion can be truncated as follows 

T B (e = E tot -e)=T B (e = E tot ). (64) 

Therefore, we understand that if Eq. i|63[l is satisfied, the environment plays the role of a heat bath because its 
temperature is almost not affected by the system energy. 

We suppose further that the interaction between the system and the environment, even if small, is able to make 
the total probability distribution microcanonical on the total energy shell at energy .Etot- We suppose also that the 
energy levels are discrete and, therefore, that E to t = E s + Ef, where s and b are discrete index's, respectively, for the 
system and the environment. Therefore, the probability Ps(E s ) for the system being at energy E s is given by 

u it? \ ^s(E fc = E tot - E B ) 

Ps(E s ) = , (65) 

"tot l^tot ) 

where Q B (Eb) is the number of states of the environment available at energy Ef,, and ^tot(Etot) the number of states 
of the total system available at energy E to t- Using Eqs. I60fl and l|63|) . one gets that 

e iS B (E b =E tot -E s ) e iS B (E b =E tot )- kTB{E E J =Etot) 
^tot(Etot) ^tot(Etot) 

Using the normalization Ps{E s ) — 1, one finally gets the well-known canonical probability distribution for the system 



p kT B (E B =E tot ) 

P S (E S ) = , (67) 



where Z = ^ s e kT B ( E B = E tat) . 



We conclude that two conditions are necessary in order to thermalize the system to a canonical probability distri- 
bution due to the contact with the environment: a large heat capacity of the environment | C vB (e)\ ^> | T ^ e ^ | and a 
microcanonical distribution on the total energy shell. 

Let us apply this result to the spin-GORM model. The density of states of the environment is 



AN 1 

n B {e) = -J--e*. (68) 
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Therefore the microcanonical temperature of the environment is 



.2 1 



and the heat capacity is 



C vB {e) = -j— r- (70) 

e + 4 

The canonical distribution of the populations of the system, at the microcanonical temperature of the environment, 
is given by 

The first condition ((63() in order to thermalize the system to a canonical probability distribution becomes 

\C vB {e)T B {e)\ = »A (72) 

e + 4 

and is depicted in Fig ED The second condition to have a microcanonical distribution on the total energy shell 
is satisfied (as we discussed in Sec. IV All when A 2 TV > 1. In this case, the populations of the system obey the 
microcanonical equilibrium value of the Markovian version of our perturbative equation 1(3 7|) . i.e., 



l-^-fF-v/l-^+f) 2 (73) 



If the two conditions are satisfied, then Eqs. 1(71(1 and 1)73(1 should be equal. The comparison between Eqs. ((71(1 and 
((7311 can be seen in Fig. ^|for different system energies. We see that the smaller the system energy is the better the 
comparison is. 

Therefore, we can conclude that under these two conditions (| 2 f | 3> A and A N > 1), the random matrices of 

e + 4 

the spin-GORM model can model an environment that behaves as a heat bath. 



C. Thermalization of the total system 



Until now, we have chosen initial conditions where the system is in the upper state with a microcanonical environ- 
ment at a given energy, like in Eqs. 1(56|) or l(57|) . We now want to consider initial conditions where the system is 
again in the upper state but where the environment is at a given canonical temperature: 



m = |i )( i| ® __ = £ -^_|i„)<i„|. (74) 

Therefore 



(a|p(0)|a)=^^^|(a|ln)| 2 . (75) 

n 

It is important to notice that there is no statistical equivalence between the canonical and the microcanonical ensembles 
in the spin-GORM model (see Appendix B). Therefore, it is interesting to ask how the probability distribution looks 
like at equilibrium after the interaction. 

One can clarify this point by plotting (a|/5(0)|a) versus energy. One uses the following energy representation 



P(e)=Y,8(E a -e)(a\p(0)\a). 

a. 



(76) 
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because 



If the total system thermalizes and reaches a canonical distribution for the total system at an effective temperature 
/3~g , one would have that 

p— A»ff£ 

P(e) = -=—, (77) 

^tot 

(«l/3(0)|a) = —z ■ (78) 

^tot 

As we shall see, it is the case if, again, A is large enough to induce "mixing" between the states X 2 N > 1. 

Indeed, one sees in Fig. ITTT a') that, for A = 0, the states of the total system corresponding to the upper level of 
the system are exponentially populated and the ones corresponding to the lower level are not. When the interaction 
is turned on and increased, one can notice that the probability distribution starts to accumulate around a mean 
effective canonical distribution. As expected, this accumulation becomes significant when X 2 N > 1 and, in this case, 
the total system can be considered as having thermalized. One can calculate the final effective temperature that the 
total system has reached after interaction. This effective temperature is depicted in Fig. I17f b) as a function of A. 
The correlation coefficient indicates whether the exponential fit of the final effective temperature is good or not. The 
effective temperature obeys the following law: f3 e f / = ± + x i ■ We also show in Fig. I17f c) the comparison between the 
time-averaged value of (a z ) and his canonical average computed with the effective temperature. One sees that, when 
A 2 A > 1, both coincide. 

This thermalization is not statistical in the sense of the equivalence between the ensembles. It is an intrinsic 
thcrmalization due to the complexity of the interaction between the states. This thermalization appears at a critical 
value of the coupling parameter when the interaction term becomes of the order of or larger than the mean level 
spacing of the total system. 



VI. CONCLUSIONS 



In this paper, we have studied a system made of two parts: a two-level system interacting in a nondiagonal way 
with a complex environment modeled by Gaussian orthogonal random matrices. 

Wc began our study by analyzing the spectral properties of this model. We investigated the spectrum on a large 
energy scale with the averaged smooth density of states and on a finer energy scale with the eigenvalue diagrams, the 
shape of the eigenstates (SOE), the spacing distribution, and the asymptotic transition probability kernel (ATPK). 
We found a global repulsion as well as avoided crossings between the eigenvalues when the coupling parameter A was 
increased. We also showed the existence of three regimes (easy to distinguish in the SOE) that are important to 
describe the different qualitative behaviors of the model: the localized regime when the interaction between the levels 
is weaker then the mean level spacing X 2 N < 1 (giving rise to very narrow SOE), the Lorentzian regime when the 
interaction between the levels becomes larger then the mean level spacing A 2 A > 1 (giving rise to a Lorentzian SOE), 
and the delocalized regime for very large A (giving rise to SOE spread over the whole spectrum). 

After the spectral study, we started the study of the dynamics of the system populations. We defined different 
domains in the parameter space (see Fig. ^) and related them to the different relaxation behaviors of the system 
population induced by the interaction with the environment. For each of these domains we tested the validity of 
approximated population evolution equations. In the strong coupling limit, we identified a population relaxation 
regime described by a Bessel function: ~ 2 that scales in time according to Xt and reaches an equilibrium 

distribution corresponding to the same probability of being in the upper and lower states of the system. In the small 
coupling limit and for small system energy, we obtained a Pauli-type equation (|36() - (|38|l describing an exponential 
relaxation of the system population that scales in time according to X 2 t and that reaches an equilibrium distribution 
depending on the system energy. The equilibrium value corresponds to a microcanonical probability distribution of 
being in a nonperturbed state of the total system inside the total energy shell. Finally, we showed the necessity of 
taking into account the non-Markovian effects in the dynamics (which are important when the system energy becomes 
non-negligible in front of the environment energy) using a non-Markovian perturbative equation i|33[) and (|34|l derived 
by the authors in Ref. |17| . This equation is perturbative and therefore only valid for small coupling parameters. This 
equation describes the highly non-Markovian dynamics of the population (made of small and fast system oscillations 
damped on a time scale corresponding to the environment time scale) when the system energy becomes large in front 
of the environment energy. This equation also reduces to the Pauli-type equation in the opposite situation, when the 
system energy is small compared to the environment energy. The validity of these approximated equations depend on 
the parameter domain considered and are summarized in Fig. 1121 An important result concerning the small coupling 
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limit is that there exist lower and upper bounds on the coupling parameter values for which the pertubative equation 
holds. The lower bound depends on the spacing between the states of the total systems and therefore on the size 
TV of the random matrices modeling the environment. We showed that this lower bound is related to the transition 
between the localized and the Lorentzian regimes in the SOE and that this bound disappears when TV — > oo. 

Another important result concerns the equilibrium values of the spin-GORM model. We showed that they are very 
well reproduced by the time averaged quantities. Moreover, we showed that the spin-GORM model can, under two 
conditions, describe the thermalization of the system to a canonical energy probability distribution corresponding 
to the environment microcanonical temperature. The two conditions for this thermalization are: a microcanonical 
probability distribution on the energy shell of the total system (that occurs when A 2 N > 1) and a large heat capacity 
of the environment compared to the ratio between the characteristic system energy and the environment temperature. 

Finally, we showed that the spin-GORM model can undergo an intrinsic thermalization due to the complex interac- 
tion between the states, and reach an overall thermal canonical distribution. This thermalization again occurs when 
the coupling parameter is large enough (i.e., larger than the mean level spacing of the total system) to "mix" the 
levels. 
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APPENDIX A: Gaussian orthogonal random matrices (GORM) 



A Gaussian orthogonal random matrix (GORM) Y is characterized by M, the size of the matrix, and by the 

parameter a^>, which enters the Gaussian probability distribution P(Y) — Ce~~J~ Tr ( y2 - 1 of the whole matrix. The 
statistical properties of a GORM are preserved under orthogonal transformations. Because the matrix is orthogonal, 
each nondiagonal element Yij is equal to its transposed Yji. The M ( J ^+ 1 ) independent matrix elements of Y are 
Gaussian random numbers of mean zero. The standard deviation of the nondiagonal matrix elements <J^ D and the 
standard deviation of the diagonal matrix element are related to Oy by 

a% = V2a^ D = J— . (79) 

The density of states of the GORM Y is defined by 

d{E)=Y,S(E-E i ), (80) 

i=l 

and the smoothed density of states by 



E+4 
E- 



d(E) = Urn - / ,1(E) .IE. (SI) 



where e is a small energy interval which is large enough to contain many states in order for d{E) to be smooth. The 
averaged smoothed density of states is an ensemble average of y realizations of the GORM. Such an ensemble is called 
the Gaussian orthogonal ensemble (GOE). It is well known |l4. Il5l ITfil| that the ensemble averaged smoothed density 
of states (d(E)) x obey the Wigner semicircle law in the limit y — *■ oo: 




(d(f7))oo - -^J E? if \E\< 



I2M 



Y 



if \E\ > J—. (82) 

V a Y 

'8M 



The domain of energy where the eigenvalues are distributed (i.e., the width of the semi-circle) is VY — y ^t"- Notice 

that when M — > oo, d(E) — * (<i(£ , )) 00 , and therefore d(E) follows the semi-circle law. The following notation is used 
in the present paper: n w (E) = (d(E)) OQ . 
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APPENDIX B: Equivalence between ensembles for the environment 



Let us consider a system interacting with its environment. The environment is in a canonical distribution at 
1 

k[j- 



temperature T can = . The evolution of a system observable is given by 



P -/3H B „ 

(A s f(t) = Tr pg^—^e^Ase-^K (83) 

On the other hand, if the environment is in a microcanonical ensemble, the evolution of the system observable is given 
by 

(A s ) £ {t) = Tr p s (0) S( - e zJM e iH M t^ se -iH tot t_ (84) 
n(H B ) 

One therefore sees that 

(AsAt) = J den{e) e -^-(A s y(t). (85) 

The statistical equivalence between the canonical and microcanonical ensembles, (As)^(t) = (As) e (t), thus occurs 
when 

n(e) e -^*5(e-e'). (86) 



One understands that this equivalence is qualitatively satisfied when n(e) is an incresing function of e. In this case, 
is a sharply peaked function. Ir 

derivative of its logarithm In n(e)^^- 



n(e)^z^- is a sharply peaked function. In order to find the maximum of n(e)^^-, we require the vanishing of the 

1 = 0. We find that 



^-S'(emax) — Tf, J- \ ~ Tf, ! (87) 

1 micro \€ max ) J- can 

where the microcanonical entropy is given by 

5(e) = k\nn(e)Se, (88) 

where <5e is a small energy shell containing many levels. This shows that if n{e)^^- is a sharply peaked function 
around e ma x, the canonical average at temperature T can is equivalent to the microcanonical average at the energy 
£max corresponding to the microcanonical temperature T m i cro (e) = T can . 

For the spin-GORM model, there is no equivalence between the canonical and the microcanonical ensembles. It is 
due to the fact that the semicircular energy distribution is not an increasing function of the energy. In this sense, this 
shows that the semicircular energy distribution does not describe a usual environment. 

To complete our reasoning, we notice that the maximum of n(e)%?^- is given by 



so that e max = and e ma x = 
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APPENDIX C: Perturbation theory 



There is no analytical way of getting a general form of the eigenvalues E a of the total system, but the three terms 
in Eq. JSJ) have different orders of magnitude depending of the value of the parameters A and A. The system and the 
environment Hamiltonians are, respectively, of order A and 1 while the coupling term is of order A. Therefore, we 
can examine the different extreme cases that can be treated perturbatively. 

A,l > A : 

When the system and the environment Hamiltonians are larger than the interaction term in Eq. we can treat 
the interaction term in a perturbative way, taking the system and the environment Hamiltonian as reference, 

H \sb) = E° sb \sb) , (90) 

where we replaced the index n by the two indices s, b. The perturbed energy is given to the second order by 

E a = E s , b = + E* + A 2 £ J (&/ ™ 2 gA + 0(A<). (91) 

b'^b b b 

We notice that the first nonzero correction to the nonperturbed eigenstate is of order A 2 . 
A » 1, A : 

When A is large compared to A and 1 in Eq. (J3J), it is possible to consider the interaction term as the reference 
Hamiltonian and to treat H$ & n d Hb as small perturbation. Transforming (JSj by a unitary matrix acting only on 
the system degree of freedom, we get 

fltot = y °x +H B + \cr z B. (92) 
The nonperturbed reference Hamiltonian is, therefore, 

H = Xa z B. (93) 
Let E KV and \nr)) = \k) \rj) be, respectively, the eigenvalues and eigenvectors of H : 

H \nri) = \Ba z \nr]) = \E Kr) \nr)) = \KE v \nr)), (94) 

where T] = 1, -y and k = ±1. The energy of the perturbed Hamiltonian is thus given to the second order perturbation 
in i by 

T = ^ + A^ 1 ^ A 2 ^ ^ - EgZ + ° UJ ' ( } 

re', 77' ' 

1 > A, A : 

In this case, the bath Hamiltonian is large compared to the system Hamiltonian and the interaction term so that 
they both can be considered as perturbations. We get 

E a = E s , b = + E* + A 2 J2 W3§F + 0(A 2 ) + 0(A 2 ). (96) 

b'=ib b b ' 

A > 1,A : 
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We now suppose that the system Hamiltonian, taken as reference, is large compared to the environment Hamiltonian 
and the interaction term, so that these last two terms can be considered as perturbations. We then get 

E a =E Stb = |l s + i? ; f + s ^^|(6'|B|6)| 2 + 0(l) + 0(A 2 ). (97) 

b'^b 

Two more situations, 1,A > A and A,A> 1 , could be considered but cannot be treated perturbativcly because 
no reference basis exists in which Hb and B are simultaneously diagonal. 
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FIG. 1: Representation of the three different domains in the space of the reduced parameters A and A of the model for a fixed 
number N of states. 
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FIG. 2: Smoothed densities of states n w (e - 
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of states of the nonperturbed spectrum is obtained by the sum of them (see Eq. 1161 '). 



(e+ v) for different values of A. The total smoothed averaged density 
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FIG. 3: Total smoothed averaged density of states obtained numerically for different values of A and A. (a) corresponds to 
A = 0.01, (b) to A = 0.5, and (c) and (d) to A = 5. For all of them N = 500 and x = 50. Notice that the A = and the 
A = 0.1 curves are not distinguishable. 




FIG. 4: Different parts of the eigenvalue diagrams with N = 500, corresponding to (a) and (b) A = 0.01, (c) and (d) A = 0.5, 
and (e) A = 5. They represent the eigenvalues of the total Hamiltonian 114H as a function of A. 
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FIG. 6: Diagram of the different regimes as a function of the coupling parameter A. 
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FIG. 7: Schematic representation of the different regimes in the plane of the reduced parameter A and N. The Lorentzian 1 
regime corresponds to T < Se and the Lorentzian 2 regime to F > Se. 




FIG. 8: (a) SOE in the perturbative regime, (b) SOE in the Lorentzian regime, (c) The width and the correlation coefficient of 
the fit of the SOE by a Lorentzian. (d) Power-law dependence of the width of the Lorentzian SOE in the coupling parameter. 
In all the figures: A = 0.01, N = 500, \ = 50, and e = 0. 
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FIG. 9: (a) ATPK from the localized regime, through the Lorentzian regime, to the delocalized regime. The ATPK is 
microcanonically averaged over the different initial conditions corresponding to the levels inside the energy shell centered at 
e° = with width 5e° = 0.05. "Pin" denotes the proportion of the ATPK that stays inside the initial energy shell after an 
infinite time, (b) and (c) ATPK for a single level as an initial condition, without any average. All the figures are obtained for 
very small system energy A = 0.01 with no random matrix ensemble average \ = 1 and for N = 500. 




FIG. 10: Time evolution ol the z component of the spin: comparison between the Markovian (M) and non-Markovian (NM) 
versions of the perturbative equation for different values of the coupling parameter. In all the cases A = 0.01 and e = 0. 
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FIG. 11: Time evolution of the z component of the spin, (a), (b), and (c) Comparison between the exact von Neumann equation 
and the non-Markovian (NM) version of the perturbative equation for different values of the coupling parameter, (a) and (b) 
correspond to A = 0.1 and (c) to A = 0.5 and A = 0.1. (d) Comparison between the exact von Neumann equation and the 
Bessel strong coupling result given by Eq. Q53p f° r A = 0.01, N — 500, and \ = 50. In all the figures \ = 50. 
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FIG. 12: Table giving the validity of the approximated equation and of the time scaling for the different parameter domains 
of the moodel. The parameter domains are defined in Fig. Q "M pert, eq." means the Markovian verion of the perturbative 
equation 13611 . "NM pert, eq." refers to the Markovian version of the perturbative equation I33L and "Bessel" refers to the 
equation I53H obtained in the very strong coupling regime. 
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FIG. 13: In all the figures A = 0.1, N = 500, e = 0, and 8e = 0.05. (a) Variance between individual trajectories and the averaged 
one (x = 100) as a function of time for different values of X 2 N. (b) Power-law dependence between the equilibrium value of this 
variance and A 2 iV. (c)-(f) Individual trajectories of the ensemble (dashed lines) and the ensemble averaged trajectory (solid 
line). In (c) X 2 N = 0.1, in (d) A 2 iV = 1, in (e) A 2 iV = 10, and in (f) X 2 N = 100. 
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FIG. 14: Comparison between between the equilibrium value of (a z )oo given by the Pauli equation and the exact values given by 
the time averaged value for different values of A. These equilibrium values are depicted as a function of the initial microcanonical 
energy e of the environment. Se — 0.05 in all figures. The parameter values are (a) A = 0.01 and N = 500; (b) A = 0.5 and 
N = 500; (c) A = 0.01 and N = 200 — 2000; (d) A = 0.5 and N = 200 - 2000. 




e 

FIG. 15: Representation of the region where the thermalization condition — — f- > A holds. 
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FIG. 16: Comparison between Eqs. 1711 and 1731 for different values of A. The narrow dotted lines are plotted to show the 
undetermination around the environment energy. The parameter values are (a) A = 0.01; (b) A = 0.05; (c) A = 0.1. "mic" 
means microcanonic and "can" means canonic. 
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FIG. 17: In all the figures A = 0.01, N = 500, and (3 = 2. (a) Probability P(e) = Dist(e) of being in an eigenstate of the total 
system at equilibrium starting from the initial condition p(0) = e g g B with i = 0.5. (b) Effective temperature of the 

equilibrium probability distribution obtained by fitting a canonical distribution to the data of (a), (c) Comparison between the 
exact equilibrium population value and the canonical one computed with the effective temperature of plot (b). 



